BACKGROUND

When describing these time series, we have used words such as “trend” and “seasonal” which need to be defined more carefully.

Trend A trend exists when there is a long-term increase or decrease in the data. It does not have to be linear. Sometimes we will refer to a trend as “changing direction”, when it might go from an increasing trend to a decreasing trend, like this figure;

SIC Trend

Seasonal A seasonal pattern occurs when a time series is affected by seasonal factors such as the time of the year or the day of the week. Seasonality is always of a fixed and known period.

All this analysis are do to SubArea 48.1 in the Southern Ocean.

OBJECTIVE

To condition an operational model in krill, the chosen environmental component must have an associated error (variance) to be incorporated, which gives more reliability to the projection of the population variables.

library(here)
library(tidyverse)
library(tsibble)
library(zoo)
library(feasts)
library(kableExtra)
library(CCAMLRGIS)

Read data. This data is from another analysis with raster layer from satellite sources.

load("DataEnvKrill2.RData")
# ls()
#cargo objeto

dtats <- get("dataenvf2")

Time Series SIC

Data wrangling to clean it, becauso have different structure and NA, wich is not supported.

dtats1 <- dtats %>% 
  select(1,2,4,6) %>% 
  mutate(ANOn =as.numeric(ANO)) %>% 
  drop_na(ANOn)

Change dtats1 from data.frameto ts object.

ts_obj <- ts(dtats1$meansic, 
             start = c(dtats1$ANOn[1], 1), 
             end = c(dtats1$ANOn[length(dtats1$ANOn)], 12), 
             frequency = 12)

ts_tbobj <- as_tsibble(ts_obj)

A subseries plot to show difference.

ts_tbobj |>
  gg_subseries(value) +
  labs(
    y = "SIC",
    title = "SIC SubArea 48.1"
  )+
  theme_bw()

### Autocorrelacion

An scaterplot of correlation simple bettwen lags. We select 60 lags, to get entire series of SIC.

ts_tbobj %>% 
  gg_lag(value, 
         geom="point",
         lags = 1:12)

ts_tbobj |> 
  ACF(value, lag_max = 60) %>% 
  autoplot() +
  labs(title = "Autocorrelacion de SIC in 48.1")+
  theme_bw()

We then describe a trend correlation and also a seasonal one, which is very characteristic of environmental data.

Descomposicion

In this chapter, we consider the most common methods for extracting these components from a time series. Often this is done to help improve understanding of the time series, but it can also be used to improve forecast accuracy.

When decomposing a time series, it is sometimes helpful to first transform or adjust the series in order to make the decomposition (and later analysis) as simple as possible. So we will begin by discussing transformations and adjustments.

dcmp <- ts_tbobj|>
  model(stl = STL(value))

otra manera de edescomponer y forecast

descomsic <- decompose(ts_obj)
#autoplot(descomsic)
acf(ts_obj)

pacf(ts_obj)

dy <-diff(ts_obj)
print(modeloarsic <- forecast::auto.arima(ts_obj,
                                       d=1,
                                       D=1,
                                       stepwise = F,
                                       approximation = F,
                                       trace=F))

## Series: ts_obj 
## ARIMA(4,1,0)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4     sma1
##       -0.5222  -0.4178  -0.1507  -0.1104  -0.8661
## s.e.   0.0467   0.0491   0.0495   0.0444   0.0360
## 
## sigma^2 = 2288:  log likelihood = -2728.67
## AIC=5469.34   AICc=5469.5   BIC=5494.8

forecast::checkresiduals(modeloarsic)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(4,1,0)(0,1,1)[12]
## Q* = 172.9, df = 19, p-value < 2.2e-16
## 
## Model df: 5.   Total lags used: 24

fctsic <- forecast::forecast(modeloarsic, h=10, level=95)
autoplot(fctsic)+
  theme_bw()

components(dcmp) |>
  as_tsibble() |>
  autoplot(value, colour="blue") +
  geom_line(aes(y=trend), colour = "#D55E00") +
  labs(
    y = "Sea Ice Concentration",
    title = "Serie de tiempo descompuesta de Sea Ice Cover en 48.1"
  ) +
  theme_bw()

Ahora lo descompongo en los demas componentes

components(dcmp) |> 
  autoplot(colour="blue") +
  theme_bw()

Model ARIMA

train <- dtats$meansic[1:40]
test <- dtats$meansic[40:43]

model1 <- arima(train, order =c(1,2,1))

predo <- forecast::forecast(model1, h=10)
plot(predo)

Moving Average.

mov_avg <- zoo::rollmean(dtats1$meansic, k=3, fill=NA)


prediction <- tail(mov_avg, n=2)

ggplot() + 
  geom_line(aes(x = as.double(dtats1$ANO), y = dtats1$meansic), 
            data = dtats, color = "blue") +
  geom_line(aes(x = as.double(dtats1$ANO), y = mov_avg), 
            data =dtats, color = "red") +
  geom_hline(yintercept = prediction, color = "green") +
  xlab("Fecha") +
  ylab("Datos") +
  ggtitle("Datos y media móvil") +
  theme_bw()

Time Series TSM

Now, we use meantsm variable

Change dtats1 from data.frameto ts object.

dtats2 <- dtats1 %>% 
  filter(ANOn>1990) %>% 
  mutate(meantsm2=meantsm-270)
ts_tsm <- ts(dtats2$meantsm2, 
             start = c(dtats2$ANOn[1], 1),
             end = c(dtats2$ANOn[length(dtats2$ANOn)], 12),
             frequency = 12)

ts_tsm2 <- as_tsibble(ts_tsm)

A subseries plot to show difference.

ts_tsm2 |>
  gg_subseries(value) +
  labs(
    y = "TSM",
    title = "TSM SubArea 48.1"
  )+
  theme_bw()

### Autocorrelacion

An scaterplot of correlation simple bettwen lags. We select 60 lags, to get entire series of SIC.

ts_tsm2 %>% 
  gg_lag(value, 
         geom="point",
         lags = 1:12)

ts_tsm2 |> 
  ACF(value, lag_max = 60) %>% 
  autoplot() +
  labs(title = "Autocorrelacion de TSM in 48.1")+
  theme_bw()

We then describe a trend correlation and also a seasonal one, which is very characteristic of environmental data.

Descomposicion

dcmptsm <- ts_tsm2|>
  model(stl = STL(value))

components(dcmptsm) |>
  as_tsibble() |>
  autoplot(value, colour="red") +
  geom_line(aes(y=trend), colour = "black") +
  labs(
    y = "TSM",
    title = "Serie de tiempo descompuesta de TSM Cover en 48.1"
  ) +
  theme_bw()

Ahora lo descompongo en los demas componentes

components(dcmptsm) |> 
  autoplot(colour="red") +
  theme_bw()

otra manera de descomponer y forecast

descom <- decompose(ts_tsm)
autoplot(descom)

acf(ts_tsm)

pacf(ts_tsm)

dy <-diff(ts_tsm)
print(modeloar <- forecast::auto.arima(ts_tsm,
                                       d=1,
                                       D=1,
                                       stepwise = F,
                                       approximation = F,
                                       trace=F))

## Series: ts_tsm 
## ARIMA(3,1,0)(2,1,0)[12] 
## 
## Coefficients:
##           ar1      ar2     ar3     sar1     sar2
##       -0.3181  -0.1663  0.1052  -0.6942  -0.2729
## s.e.   0.0575   0.0616  0.0631   0.0576   0.0641
## 
## sigma^2 = 0.04078:  log likelihood = 66.47
## AIC=-120.93   AICc=-120.7   BIC=-97.44

forecast::checkresiduals(modeloar)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,0)(2,1,0)[12]
## Q* = 302.61, df = 19, p-value < 2.2e-16
## 
## Model df: 5.   Total lags used: 24

fct <- forecast::forecast(modeloar, h=10, level=95)
autoplot(fct)+
  theme_bw()

Model ARIMA

traintsm <- dtats2$meantsm2[1:32]
testtsm <- dtats$meansic[33:40]

modeltsm <- arima(traintsm, order =c(1,2,1))

predotsm <- forecast::forecast(modeltsm, h=10)
plot(predotsm)

Moving Average.

mov_avgtsm <- zoo::rollmean(dtats2$meantsm2, k=3, fill=NA)


predictiontsm <- tail(mov_avgtsm, n=5)

ggplot() + 
  geom_line(aes(x = dtats2$ANOn, y = dtats2$meantsm2), 
            data = dtats, color = "blue") +
  geom_line(aes(x = dtats2$ANOn, y = mov_avgtsm), 
            data =dtats, color = "red") +
  geom_hline(yintercept = predictiontsm[1], color = "green") +
  xlab("Fecha") +
  ylab("Datos") +
  ggtitle("Datos y media móvil") +
  theme_bw()

Time Series Analisis Env Data Krill

Mauricio Mardones

02 March, 2023

BACKGROUND

OBJECTIVE

Time Series SIC

Descomposicion

Model ARIMA

Moving Average.

Time Series TSM

Descomposicion

Model ARIMA

Moving Average.